Effects of Neighborhood Environmental Circumstances on Thalamocortical Connectivity

Prepare data

Glasser parcel names

#Glasser regions with corresponding labels and tract names
glasser.regions <- read.csv("/cbica/projects/thalamocortical_development/Maps/parcellations/surface/glasser360_regionlist.csv") #parcel name, label name 
glasser.regions$tract <- paste0("thalamus_", glasser.regions$orig_parcelname) %>% gsub("_ROI", "_autotrack", .) %>% gsub("-", "_", .) 

S-A axis

S.A.axis.cifti <- read_cifti("/cbica/projects/thalamocortical_development//Maps/S-A_ArchetypalAxis/FSLRVertex/SensorimotorAssociation_Axis_parcellated/SensorimotorAssociation.Axis.Glasser360.pscalar.nii") #S-A_ArchetypalAxis repo
S.A.axis <- data.frame(SA.axis = rank(S.A.axis.cifti$data), orig_parcelname = names(S.A.axis.cifti$Parcel))
S.A.axis <- merge(S.A.axis, glasser.regions, by = "orig_parcelname")
S.A.axis$SA.axis.bin <- as.numeric(cut2(S.A.axis$SA.axis, g = 5)) #bins for quintile enrichment analysis
SA.colorvector <- c("#FFC228", "#FFD16A", "#e9dcf2", "#BE93C4", "#6F1382") #colors to use for plotting for each of the 5 quintiles

S.A.axis$SA.axis.quartile <- as.numeric(cut2(S.A.axis$SA.axis, g = 4)) #quartiles for sensitivity enrichment
S.A.axis$SA.axis.decile <- as.numeric(cut2(S.A.axis$SA.axis, g = 10)) #deciles for sensitivity enrichment

Environment statistics

setwd("/cbica/projects/thalamocortical_development/thalamocortical_results/environment_results/") #gam output dir, from gam_functions/fit_envGams.R

#list all .csv files in the environment results directory
files <- list.files(getwd(), pattern = '.csv', ignore.case = T, full.names = F) 

for(i in 1:length(files)){
  filename <- files[i]
  Rname <- gsub('.{4}$', '', filename)
  x <- read.csv(filename, header = TRUE)
  assign(Rname, x) 
}
rm(x)

Spatial permutation (spin) test parcel rotation matrix

perm.id.full <- readRDS("/cbica/projects/thalamocortical_development/software/rotate_parcellation/glasser_sphericalrotations_N10000.rds") #10,000 spatial null spins
spin.parcels <- glasser.regions #order of complete set of glasser parcels for spinning

Dataset-specific tract lists

#generated by tract_measures/tractlists/thalamocortical_tractlists.R
tractlist.pnc <- read.csv("/cbica/projects/thalamocortical_development/thalamocortical_structuralconnectivity/individual/PNC/PNC_tractlist_primary.csv") 
tractlist.hcpd <- read.csv("/cbica/projects/thalamocortical_development/thalamocortical_structuralconnectivity/individual/HCPD/HCPD_tractlist_primary.csv") 

Associations Between Neighborhood and Household SES and Thalamocortical Connectivity

Number of significant environmental effects

##Function to calculate the number and percent of thalamocortical connections showing a significant association with an environmental characteristic
sig.effects <- function(df, env.measure){
  enveffects.df <- df
  enveffects.df$significant <- p.adjust(enveffects.df$GAM.cov.pvalue, method = c("fdr")) < 0.05
  sigeffect.totaln <- sum(enveffects.df$significant)
  sigeffect.percent <- round(sigeffect.totaln/length(enveffects.df$significant)*100, 2)
  print(sprintf("%s thalamocortical connections (%s percent) show a significant association between thalamocortical FA and the %s", sigeffect.totaln, sigeffect.percent, env.measure))
}

PNC

Significant neighborhood-level environment effects (neighborhood SES)

sig.effects(df = envSES_maineffects_FA_glasser_pnc, env.measure = "neighborhood environment (envSES factor)")
## [1] "128 thalamocortical connections (57.4 percent) show a significant association between thalamocortical FA and the neighborhood environment (envSES factor)"

controlling for parental education

sig.effects(df = envSES_maineffects_educontrolled_FA_glasser_pnc, env.measure = "neighborhood environment (envSES factor)")
## [1] "125 thalamocortical connections (56.05 percent) show a significant association between thalamocortical FA and the neighborhood environment (envSES factor)"

Percent of positive significant neighborhood-level environment effects

(envSES_maineffects_FA_glasser_pnc %>% mutate(significant = p.adjust(GAM.cov.pvalue, method = c("fdr")) < 0.05) %>% filter(significant == TRUE) %>% filter (GAM.cov.tvalue > 0) %>% nrow() / 128) * 100
## [1] 92.96875

Brain map of neighborhood-level environment effects

envSES_maineffects_FA_glasser_pnc <- envSES_maineffects_FA_glasser_pnc %>% 
  mutate(significant = p.adjust(envSES_maineffects_FA_glasser_pnc$GAM.cov.pvalue, method = c("fdr")) < 0.05)

sigeffects.df <- envSES_maineffects_FA_glasser_pnc %>% filter(significant == TRUE)

sigeffects.positivet.df <- sigeffects.df %>% 
  filter(GAM.cov.tvalue > 0) %>% 
  left_join(glasser.regions, ., by = "tract") %>%
  select(label, tract, GAM.cov.tvalue)

sigeffects.positivet.df$tractlist <- sigeffects.positivet.df$tract %in% tractlist.pnc$tract

ggplot() + 
  geom_brain(data = sigeffects.positivet.df, atlas = glasser, mapping = aes(fill = GAM.cov.tvalue, colour=I("black"),  size=I(.05))) +  
  theme_void()  + 
  scale_fill_gradientn(colours = c(alpha("#323280", 0.15), "#672975"), limits = c(2, 5.25), oob = squish, na.value = "white") + 
   new_scale_fill() + 
  geom_brain(data = sigeffects.positivet.df, atlas = glasser, mapping = aes(fill = tractlist, colour=I("black"), size=I(.05))) + 
  scale_fill_manual(values = c("gray89", alpha("white", 0)), na.value = "white") +
   theme(
   legend.text = element_text(size = 5, family = "Arial", color = c("black")),
   legend.key.size = unit(1, 'mm'),
   legend.title = element_text(size = 5, family = "Arial", color = c("black")))

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure7/envSES_corticalmap.pdf", device = "pdf", dpi = 500,  width = 4.35, height = 4)

Significant household-level environment effects (average parental years of education)

sig.effects(df = householdSES_maineffects_FA_glasser_pnc, env.measure = "household environment (average parental years of education)")
## [1] "16 thalamocortical connections (7.17 percent) show a significant association between thalamocortical FA and the household environment (average parental years of education)"

controlling for envSES

sig.effects(df = householdSES_maineffects_envSEScontrolled_FA_glasser_pnc, env.measure = "household environment (average parental years of education)")
## [1] "0 thalamocortical connections (0 percent) show a significant association between thalamocortical FA and the household environment (average parental years of education)"

HCPD

Significant household-level environment effects (average parental years of education)

sig.effects(df = householdSES_maineffects_linear_FA_glasser_hcpd, env.measure = "household environment (average parental years of education)")
## [1] "0 thalamocortical connections (0 percent) show a significant association between thalamocortical FA and the household environment (average parental years of education)"

Significant household-level environment effects (state-adjusted income-to-needs)

sig.effects(df = INRstate_maineffects_FA_glasser_hcpd, env.measure = "household environment (INR)")
## [1] "0 thalamocortical connections (0 percent) show a significant association between thalamocortical FA and the household environment (INR)"

Thalamocortical Developmental Trajectories Stratified by Neighborhood Environment Factor Score

Trajectories in S-A axis quintiles

#Function to plot predicted developmental trajectories for 10th and 90th percentile envSES factor scores for bins of the S-A axis
plot_trajectories_byenv_SAaxis <- function(quintile, break.points, y.limits){
  
  # Compute average developmental trajectories for low and high factors scores within each S-A axis quintile
  df <- devtrajectories_byenvSES_FA_glasser_pnc %>% merge(.,  S.A.axis, by = c("tract"))
  df$age <- round(df$age, 3)
  df$envSES <- as.factor(df$envSES)
  df <- envSES_maineffects_FA_glasser_pnc %>% select(tract, significant) %>% merge(., df, by = "tract")
  df <- df %>% filter(significant == TRUE)
  
  SAaxismooths.by.envSES <- df %>% group_by(envSES, SA.axis.bin, age) %>% 
                         do(est.mean = mean(.$fitted)) %>% 
                         unnest(cols = c(est.mean)) 

  
  # Plot the trajectories for the requested quintile
  SAtrajectories.byenv.plot <- SAaxismooths.by.envSES %>% filter(SA.axis.bin == quintile) %>%
  ggplot(aes(x = age, y = est.mean, group = interaction(envSES, SA.axis.bin), linetype = envSES)) +
  geom_line(color = SA.colorvector[quintile], linewidth = 0.4) +
  theme_classic() +
  labs(x = "\nAge (years)", y = "\n") +
  scale_y_continuous(breaks = break.points, limits = y.limits) +
  theme(
        legend.position = "none", 
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.text = element_text(size = 5, family = "Arial", color = c("black")),
        axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
        axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
        axis.line = element_line(linewidth = .2), 
        axis.ticks = element_line(linewidth = .2))
  
  return(SAtrajectories.byenv.plot)  
}
left_join(tractlist.pnc, S.A.axis, by = "tract") %>% filter(SA.axis.bin == 1) %>%
  mutate(SA.axis.bin = as.factor(SA.axis.bin)) %>%
  ggplot() + 
  geom_brain(atlas = glasser, mapping = aes(fill = SA.axis.bin, colour=I("black"),  size=I(.05))) +
  theme_void()  + 
  scale_fill_manual(values = SA.colorvector[1], na.value = "white")
## merging atlas and data by 'label'

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure7/SAquintile1.pdf", device = "pdf", dpi = 500, width = 2.3 , height = 3)
## merging atlas and data by 'label'
plot_trajectories_byenv_SAaxis(quintile = 1, break.points = c(0.41, 0.42, 0.43, 0.44), y.limits = c(0.405, 0.44))

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure7/Devtrajectory_byenvSES_quintile1.pdf", device = "pdf", dpi = 500, width = 1.27 , height = 1)
left_join(tractlist.pnc, S.A.axis, by = "tract") %>% filter(SA.axis.bin == 2) %>%
  mutate(SA.axis.bin = as.factor(SA.axis.bin)) %>%
  ggplot() + 
  geom_brain(atlas = glasser, mapping = aes(fill = SA.axis.bin, colour=I("black"),  size=I(.05))) +
  theme_void()  + 
  scale_fill_manual(values = SA.colorvector[2], na.value = "white")
## merging atlas and data by 'label'

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure7/SAquintile2.pdf", device = "pdf", dpi = 500, width = 2.3 , height = 3)
## merging atlas and data by 'label'
plot_trajectories_byenv_SAaxis(quintile = 2, break.points = c(0.37, 0.38, 0.39, 0.40), y.limits = c(0.365, 0.40))

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure7/Devtrajectory_byenvSES_quintile2.pdf", device = "pdf", dpi = 500, width = 1.27 , height = 1)
left_join(tractlist.pnc, S.A.axis, by = "tract") %>% filter(SA.axis.bin == 3) %>%
  mutate(SA.axis.bin = as.factor(SA.axis.bin)) %>%
  ggplot() + 
  geom_brain(atlas = glasser, mapping = aes(fill = SA.axis.bin, colour=I("black"),  size=I(.05))) +
  theme_void()  + 
  scale_fill_manual(values = "#d4bfe3", na.value = "white")
## merging atlas and data by 'label'

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure7/SAquintile3.pdf", device = "pdf", dpi = 500, width = 2.3 , height = 3)
## merging atlas and data by 'label'
plot_trajectories_byenv_SAaxis(quintile = 3, break.points = c(0.33, 0.34, 0.35, 0.36), y.limits = c(0.325, 0.36))

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure7/Devtrajectory_byenvSES_quintile3.pdf", device = "pdf", dpi = 500, width = 1.27 , height = 1)
left_join(tractlist.pnc, S.A.axis, by = "tract") %>% filter(SA.axis.bin == 4) %>%
  mutate(SA.axis.bin = as.factor(SA.axis.bin)) %>%
  ggplot() + 
  geom_brain(atlas = glasser, mapping = aes(fill = SA.axis.bin, colour=I("black"),  size=I(.05))) +
  theme_void()  + 
  scale_fill_manual(values = SA.colorvector[4], na.value = "white")
## merging atlas and data by 'label'
## Warning: Some data not merged properly. Check for naming errors in data:
##   atlas type  hemi  region side  label           geometry tract                 
##   <chr> <chr> <chr> <chr>  <chr> <chr>     <MULTIPOLYGON> <chr>                 
## 1 <NA>  <NA>  <NA>  <NA>   <NA>  lh_L_10pp          EMPTY thalamus_L_10pp_autot…
## # ℹ 5 more variables: orig_parcelname <chr>, SA.axis <dbl>, SA.axis.bin <fct>,
## #   SA.axis.quartile <dbl>, SA.axis.decile <dbl>
## Warning: Some data not merged. Check for spelling mistakes in:
##          label                     tract orig_parcelname SA.axis SA.axis.bin
## 396 lh_L_10pp thalamus_L_10pp_autotrack      L_10pp_ROI     282           4
##     SA.axis.quartile SA.axis.decile
## 396                4              8

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure7/SAquintile4.pdf", device = "pdf", dpi = 500, width = 2.3 , height = 3)
## merging atlas and data by 'label'
## Warning: Some data not merged properly. Check for naming errors in data:
##   atlas type  hemi  region side  label           geometry tract                 
##   <chr> <chr> <chr> <chr>  <chr> <chr>     <MULTIPOLYGON> <chr>                 
## 1 <NA>  <NA>  <NA>  <NA>   <NA>  lh_L_10pp          EMPTY thalamus_L_10pp_autot…
## # ℹ 5 more variables: orig_parcelname <chr>, SA.axis <dbl>, SA.axis.bin <fct>,
## #   SA.axis.quartile <dbl>, SA.axis.decile <dbl>

## Warning: Some data not merged. Check for spelling mistakes in:
##          label                     tract orig_parcelname SA.axis SA.axis.bin
## 396 lh_L_10pp thalamus_L_10pp_autotrack      L_10pp_ROI     282           4
##     SA.axis.quartile SA.axis.decile
## 396                4              8
plot_trajectories_byenv_SAaxis(quintile = 4, break.points = c(0.33, 0.34, 0.35, 0.36), y.limits = c(0.33, 0.36))

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure7/Devtrajectory_byenvSES_quintile4.pdf", device = "pdf", dpi = 500, width = 1.27 , height = 1)
left_join(tractlist.pnc, S.A.axis, by = "tract") %>% filter(SA.axis.bin == 5) %>%
  mutate(SA.axis.bin = as.factor(SA.axis.bin)) %>%
  ggplot() + 
  geom_brain(atlas = glasser, mapping = aes(fill = SA.axis.bin, colour=I("black"),  size=I(.05))) +
  theme_void()  + 
  scale_fill_manual(values = SA.colorvector[5], na.value = "white")
## merging atlas and data by 'label'

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure7/SAquintile5.pdf", device = "pdf", dpi = 500, width = 2.3 , height = 3)
## merging atlas and data by 'label'
plot_trajectories_byenv_SAaxis(quintile = 5, break.points = c(0.34, 0.35, 0.36, 0.37), y.limits = c(0.34, 0.375))

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure7/Devtrajectory_byenvSES_quintile5.pdf", device = "pdf", dpi = 500, width = 1.24 , height = 1)

Neighborhood Environment Effects Vary Across the Sensorimotor-Association Axis

S-A axis correlation

Spearman’s rho

spin.data <- left_join(spin.parcels, sigeffects.df, by = "tract") #merged data for spin tests, with all 360 parcels in the expected right hemisphere --> left hemisphere order (matching glasser.coords)
spin.data <- merge(spin.data, S.A.axis, by = c("label", "orig_parcelname", "tract"), sort =F)

cor.test(spin.data$GAM.cov.tvalue, spin.data$SA.axis, method = c("spearman"))
## 
##  Spearman's rank correlation rho
## 
## data:  spin.data$GAM.cov.tvalue and spin.data$SA.axis
## S = 244280, p-value = 0.0005847
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.3010667

Spatial permutation (spin) based p-value

perm.sphere.p(x = spin.data$SA.axis, y = spin.data$GAM.cov.tvalue, perm.id = perm.id.full, corr.type = "spearman") 
## [1] 0.0056

Correlation plot

ggplot() +
  geom_point(data = spin.data, aes(x = SA.axis, y = GAM.cov.tvalue, color = SA.axis), size = 0.5) +
  geom_smooth(data = spin.data, aes(x = SA.axis, y = GAM.cov.tvalue), method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) + 
  scale_y_continuous(limits = c(.5, 10)) +
  scale_color_gradient2(low = "goldenrod1", mid = "#ede4f5", high = "#6f1282", guide = "colorbar", aesthetics = "color", midpoint = 180) +
  theme_classic() +
  labs(x="\nS-A axis", y="Environment effect (t)") +
  theme(
  legend.position = "none", 
  axis.text = element_text(size = 5, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2)) 

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure7/EnvEffect_SArank_correlationplot.pdf", device = "pdf", dpi = 500, width = 1.61, height = 1.3)

Specificity test: cortical and thalamic axis correlations

#Glasser parcel x,y,z coordinates
load(file = "/cbica/projects/thalamocortical_development/Maps/parcellations/surface/hcp_mmp1.0.rda") #RDA file from the brainGraph package https://github.com/cwatson/brainGraph/blob/master/data/hcp_mmp1.0.rda with glasser region x, y, and z MNI coordinates. regions are in lh --> rh order
hcp_mmp1.0$x.mni[1:180]  <- hcp_mmp1.0$x.mni[1:180]*-1 #convert R-> L hemisphere coords into medial lateral coords
hcp_mmp1.0$z.mni <- hcp_mmp1.0$z.mni*-1 
hcp_mmp1.0$label <- NA
hcp_mmp1.0$label[1:180] <- glasser.regions$label[181:360]
hcp_mmp1.0$label[181:360] <- glasser.regions$label[1:180]
hcp_mmp1.0 <- hcp_mmp1.0 %>% select(-hemi)
#C-Mt files for defining the group-level core-matrix gradient
CPt.glasser.pnc <- read.csv("/cbica/projects/thalamocortical_development/thalamocortical_structuralconnectivity/individual/PNC/PNC_tractaverage_CPt_primary.csv") %>% select("tract", "label", "orig_parcelname", "mean_CPt")
CPt.glasser.hcpd <- read.csv("/cbica/projects/thalamocortical_development/thalamocortical_structuralconnectivity/individual/HCPD/HCPD_tractaverage_CPt_primary.csv")  %>% select("tract", "label", "orig_parcelname", "mean_CPt")
CPt.gradient <- merge(CPt.glasser.pnc, CPt.glasser.hcpd, by = c("label", "tract", "orig_parcelname"), all.x = TRUE, all.y = TRUE, suffixes = c("_pnc", "_hcpd"))
CPt.gradient$CPt <- CPt.gradient %>% select(mean_CPt_pnc, mean_CPt_hcpd) %>% rowMeans(na.rm = T)
CPt.gradient <- CPt.gradient %>% select(CPt, label)

hcp_mmp1.0 <- left_join(hcp_mmp1.0, CPt.gradient, by = "label")
#Function to compute the correlation between an environment map and an anatomical axis (x,y,z) as well as test the significance of the difference between the correlation with the anatomical axis versus the S-A axis
compute_axis_correlation <- function(axis){
  df <- merge(spin.data, hcp_mmp1.0, by = "label", sort = F)
  completeobs.n <- sum(!is.na(df$GAM.cov.tvalue))
  
  #S-A axis - environment correlation (as computed above) for comparison of two overlapping correlations based on dependent groups
  SA.axis.cor <-  cor(df$SA.axis, df$GAM.cov.tvalue, use = "complete.obs", method = c("spearman"))
  
  #Anatomical axis - environment correlation
  anatomical.axis.cor <- cor(df[,axis], df$GAM.cov.tvalue, use = "complete.obs", method = c("spearman")) #correlation between anatomical axis and env metric
  anatomical.axis.pvalue <- perm.sphere.p(x = df[,axis], y = df$GAM.cov.tvalue, perm.id = perm.id.full, corr.type = "spearman") #spin based p-value for the correlation
  
  #S-A axis - anatomical axis correlation
  SA.anatomical.cor <- cor(df$SA.axis, df[,axis], use = "complete.obs", method = c("spearman"))
  
  #Test for significant difference between two correlations based on dependent groups (i.e., correlations with an overlapping input)
  cocor.result <- cocor.pvalue <- cocor.dep.groups.overlap(
  r.jk = SA.axis.cor, 
  r.jh = anatomical.axis.cor, 
  r.kh = SA.anatomical.cor, 
  n = completeobs.n, alternative = "two.sided", test = "hittner2003")
  
  cocor.pvalue <- cocor.result@hittner2003[["p.value"]]
  
  comparison.results <- data.frame(axis = axis, axis.cor = anatomical.axis.cor, axis.cor.pvalue = anatomical.axis.pvalue, SA.cor = SA.axis.cor, cocor.pvalue = cocor.pvalue)
  
  return(comparison.results)
  }

Anterior-posterior axis

compute_axis_correlation(axis = "y.mni")
##    axis   axis.cor axis.cor.pvalue    SA.cor cocor.pvalue
## 1 y.mni 0.06846274         0.38685 0.3010667 0.0005890574

Medial-lateral axis

compute_axis_correlation(axis = "x.mni")
##    axis  axis.cor axis.cor.pvalue    SA.cor cocor.pvalue
## 1 x.mni 0.1841123         0.10335 0.3010667    0.3319559

Dorsal-ventral axis

compute_axis_correlation(axis = "z.mni")
##    axis  axis.cor axis.cor.pvalue    SA.cor cocor.pvalue
## 1 z.mni 0.2117744         0.07605 0.3010667    0.3909002

Core-matrix gradient

compute_axis_correlation(axis = "CPt")
##   axis   axis.cor axis.cor.pvalue    SA.cor cocor.pvalue
## 1  CPt -0.1265908         0.15965 0.3010667   4.3867e-07

Correct axis correlation ps

axis.cor.ps <- c(0.0056, 0.38685, 0.10335, 0.07605, 0.15965) #S-A, A-P, M-L, D-V, C-M
p.adjust(axis.cor.ps, method = c("fdr"))
## [1] 0.0280000 0.3868500 0.1722500 0.1722500 0.1995625

Correct correlation comparison ps

anatomical.ps <- c(0.0005890574, 0.3319559, 0.3909002, 4.3867e-07) #A-P, M-L, D-V, C-M
p.adjust(anatomical.ps, method = c("fdr"))
## [1] 1.178115e-03 3.909002e-01 3.909002e-01 1.754680e-06
#Put above results into df for plotting
axis.results <- data.frame(Axis = factor(), corr = double())
axis.results <- axis.results %>% add_row(Axis = "S-A", corr = 0.3010667)
axis.results <- axis.results %>% add_row(Axis = "A-P", corr = 0.06846274)
axis.results <- axis.results %>% add_row(Axis = "D-V", corr = 0.2117744)
axis.results <- axis.results %>% add_row(Axis = "M-L", corr = 0.1841123)
axis.results <- axis.results %>% add_row(Axis = "C-M", corr = -0.1265908)

axis.results$Axis <- factor(axis.results$Axis, ordered = TRUE, levels = c("S-A", "D-V", "M-L", "A-P", "C-M"))

ggplot(axis.results, aes(x = Axis, y = corr, fill = Axis)) + 
  geom_bar(stat='identity') + 
  labs(x="") +
  labs(y="\nEnvironment Effect Correlation") +
  theme_classic()+
  scale_fill_manual(values=c(alpha("#323280", 0.5), "#672975", "#672975", "#672975", "#672975")) +
  theme(
  legend.position = "none", 
  axis.text = element_text(size = 5, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2))

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure7/EnvEffect_anatomicalaxes_barplot.pdf", device = "pdf", dpi = 500, width = 1.58, height = 1.3)

Effect magnitude enrichment testing (S-A axis quintiles)

Mean environment effect t-value by S-A axis quintile (empirical)

#Calculate empirical enrichment statistic for each S-A axis bin
envSES_effects_SAaxisbins <- spin.data %>% 
                             group_by(SA.axis.bin) %>% 
                             do(mean.tvalue.significant = mean(.$GAM.cov.tvalue, na.rm = TRUE)) %>%
                             unnest(cols = c("mean.tvalue.significant"))
envSES_effects_SAaxisbins
## # A tibble: 5 × 2
##   SA.axis.bin mean.tvalue.significant
##         <dbl>                   <dbl>
## 1           1                    2.70
## 2           2                    3.60
## 3           3                    2.79
## 4           4                    3.91
## 5           5                    4.57

Null distribution and enrichment testing by S-A axis quintile

#Calculate null distribution of t-values by bin based on spherical spatial permutations (spins)

## inputs to spatially permute (i.e., spin) with the pre-computed spatial permutation matrix 
x <- spin.data$SA.axis.bin #cortical map 1 to spatially permute: S-A axis bins
y <- spin.data$GAM.cov.tvalue #cortical map 2 to spatially permute: significant envSES-thalamocortical FA t-values
perm.id <- perm.id.full

## spin the empirical data 
nroi = dim(perm.id)[1]  #number of regions
nperm = dim(perm.id)[2] #number of permutations
x.perm = y.perm = array(NA,dim=c(nroi,nperm)) #initialize
  for (r in 1:nperm) {
    for (i in 1:nroi) {
      x.perm[i,r] = x[perm.id[i,r]] #spinning x, spatially permuted bins
      y.perm[i,r] = y[perm.id[i,r]] #spinning y, spatially permuted t-values
    }
  }

## compute mean t-value for each S-A axis bin using spatially permuted y inputs 
### set up spun (y) and empirical (x) df
y.spatialperm <- as.data.frame(y.perm) %>% set_names(sprintf("perm%s", seq(from = 1, to = ncol(y.perm))))
y.spatialperm.empiricalx <- cbind(x, y.spatialperm) %>% as.data.frame()
colnames(y.spatialperm.empiricalx)[1] <- c("empirical.x")

### compute mean t-values for S-A axis bins for all spins
envSES_effects_SAaxisbins_permutedy <- y.spatialperm.empiricalx %>% group_by(empirical.x) %>% #group by empirical S-A axis bins
                                       dplyr::summarize(across(everything(), \(x) mean(x, na.rm = TRUE))) %>% #calculate mean t-value for permuted data
                                       select(-empirical.x) %>% t() %>% as.data.frame() %>%
                                       set_names(c("SA.bin1.tvalue","SA.bin2.tvalue", "SA.bin3.tvalue", "SA.bin4.tvalue", "SA.bin5.tvalue"))
#Function to plot the null distribution of environment effect t-values for a given S-A axis bin and test whether the true t-value is significantly smaller or greater in magnitude than the null (permutation based p)
enveffects_SAbin_enrichment <- function(bin, enrichment_type){
  
  nperm = 10000
  
  SAbin.tvalue.empirical <- envSES_effects_SAaxisbins %>% filter(SA.axis.bin == bin) %>% select(mean.tvalue.significant)
  
  #Histogram
  SA.colorvector <- c("#FFC228", "#FFD16A", "#e9dcf2", "#BE93C4", "#6F1382")
  
  null.tvalue.histogram <- envSES_effects_SAaxisbins_permutedy %>% select(sprintf("SA.bin%s.tvalue", bin)) %>%
    ggplot(aes(x = .[,1])) + geom_histogram(fill = c(SA.colorvector)[bin]) +
    theme_classic() +
    geom_vline(xintercept =  SAbin.tvalue.empirical$mean.tvalue.significant, linewidth = 0.25) +
    xlab("Environment effect (t-value)") +
    ylab("Number of nulls") +
    scale_x_continuous(breaks = c(2, 3, 4, 5), limits = c(1.25, 5.25)) +
    theme(
     legend.position = "none", 
     axis.text.x = element_text(size = 5, family = "Arial", color = c("black")),
     axis.text.y = element_blank(),
     axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
     axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
     axis.line.x = element_line(linewidth = .2), 
     axis.line.y = element_blank(),
     axis.ticks = element_line(linewidth = .2))
    
  #Permutation-based p
    if(enrichment_type == "smaller"){
      perm.p <- sum(envSES_effects_SAaxisbins_permutedy[,bin] < SAbin.tvalue.empirical$mean.tvalue.significant)/nperm
    }
  
  if(enrichment_type == "larger"){
      perm.p <- sum(envSES_effects_SAaxisbins_permutedy[,bin] > SAbin.tvalue.empirical$mean.tvalue.significant)/nperm
  }
  
  permutation.results <- list(null.tvalue.histogram, perm.p)
  return(permutation.results)
}

S-A quintile 1

enveffects_SAbin_enrichment(bin = 1, enrichment_type = "smaller")
## [[1]]

## 
## [[2]]
## [1] 0.0921
ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure7/quintile1_enveffect_spins.pdf", device = "pdf", dpi = 500, width = .95, height =.9)

S-A quintile 2

enveffects_SAbin_enrichment(bin = 2, enrichment_type = "smaller")
## [[1]]

## 
## [[2]]
## [1] 0.5832
ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure7/quintile2_enveffect_spins.pdf", device = "pdf", dpi = 500, width = .95, height =.9)

S-A quintile 3

enveffects_SAbin_enrichment(bin = 3, enrichment_type = "smaller")
## [[1]]

## 
## [[2]]
## [1] 0.0687
ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure7/quintile3_enveffect_spins.pdf", device = "pdf", dpi = 500, width = .95, height =.9)

S-A quintile 4

enveffects_SAbin_enrichment(bin = 4, enrichment_type = "larger")
## [[1]]

## 
## [[2]]
## [1] 0.1859
ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure7/quintile4_enveffect_spins.pdf", device = "pdf", dpi = 500, width = .95, height =.9)

S-A quintile 5

enveffects_SAbin_enrichment(bin = 5, enrichment_type = "larger")
## [[1]]

## 
## [[2]]
## [1] 0.025
ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure7/quintile5_enveffect_spins.pdf", device = "pdf", dpi = 500, width = .95, height =.9)

Enrichment sensitivity test: deciles

#Calculate empirical enrichment statistic for each S-A axis bin
envSES_effects_SAaxisbins <- spin.data %>% 
                             group_by(SA.axis.decile) %>% 
                             do(mean.tvalue.significant = mean(.$GAM.cov.tvalue, na.rm = TRUE)) %>%
                             unnest(cols = c("mean.tvalue.significant"))
#Calculate null distribution of t-values by bin based on spherical spatial permutations (spins)

## inputs to spatially permute (i.e., spin) with the pre-computed spatial permutation matrix 
x <- spin.data$SA.axis.decile #cortical map 1 to spatially permute: S-A axis bins
y <- spin.data$GAM.cov.tvalue #cortical map 2 to spatially permute: significant envSES-thalamocortical FA t-values
perm.id <- perm.id.full

## spin the empirical data 
nroi = dim(perm.id)[1]  #number of regions
nperm = dim(perm.id)[2] #number of permutations
x.perm = y.perm = array(NA,dim=c(nroi,nperm)) #initialize
  for (r in 1:nperm) {
    for (i in 1:nroi) {
      x.perm[i,r] = x[perm.id[i,r]] #spinning x, spatially permuted bins
      y.perm[i,r] = y[perm.id[i,r]] #spinning y, spatially permuted t-values
    }
  }

## compute mean t-value for each S-A axis bin using spatially permuted y inputs 
### set up spun (y) and empirical (x) df
y.spatialperm <- as.data.frame(y.perm) %>% set_names(sprintf("perm%s", seq(from = 1, to = ncol(y.perm))))
y.spatialperm.empiricalx <- cbind(x, y.spatialperm) %>% as.data.frame()
colnames(y.spatialperm.empiricalx)[1] <- c("empirical.x")

### compute mean t-values for S-A axis bins for all spins
envSES_effects_SAaxisbins_permutedy <- y.spatialperm.empiricalx %>% group_by(empirical.x) %>% #group by empirical S-A axis bins
                                       dplyr::summarize(across(everything(), \(x) mean(x, na.rm = TRUE))) %>% #calculate mean t-value for permuted data
                                       select(-empirical.x) %>% t() %>% as.data.frame() %>%
                                       set_names(c("SA.bin1.tvalue","SA.bin2.tvalue", "SA.bin3.tvalue", "SA.bin4.tvalue", "SA.bin5.tvalue", "SA.bin6.tvalue", "SA.bin7.tvalue", "SA.bin8.tvale", "SA.bin9.tvalue", "SA.bin10.tvalue"))
enveffects_SAbin_enrichment <- function(bin, enrichment_type){
  
  nperm = 10000
  
  SAbin.tvalue.empirical <- envSES_effects_SAaxisbins %>% filter(SA.axis.decile == bin) %>% select(mean.tvalue.significant)
  
  #Permutation-based p
    if(enrichment_type == "smaller"){
      perm.p <- sum(envSES_effects_SAaxisbins_permutedy[,bin] < SAbin.tvalue.empirical$mean.tvalue.significant)/nperm
    }
  
  if(enrichment_type == "larger"){
      perm.p <- sum(envSES_effects_SAaxisbins_permutedy[,bin] > SAbin.tvalue.empirical$mean.tvalue.significant)/nperm
  }
  
  return(perm.p)
}

S-A decile 10

enveffects_SAbin_enrichment(bin = 10, enrichment_type = "larger")
## [1] 0.0085

Enrichment sensitivity test: quartiles

#Calculate empirical enrichment statistic for each S-A axis bin
envSES_effects_SAaxisbins <- spin.data %>% 
                             group_by(SA.axis.quartile) %>% 
                             do(mean.tvalue.significant = mean(.$GAM.cov.tvalue, na.rm = TRUE)) %>%
                             unnest(cols = c("mean.tvalue.significant"))
#Calculate null distribution of t-values by bin based on spherical spatial permutations (spins)

## inputs to spatially permute (i.e., spin) with the pre-computed spatial permutation matrix 
x <- spin.data$SA.axis.quartile #cortical map 1 to spatially permute: S-A axis bins
y <- spin.data$GAM.cov.tvalue #cortical map 2 to spatially permute: significant envSES-thalamocortical FA t-values
perm.id <- perm.id.full

## spin the empirical data 
nroi = dim(perm.id)[1]  #number of regions
nperm = dim(perm.id)[2] #number of permutations
x.perm = y.perm = array(NA,dim=c(nroi,nperm)) #initialize
  for (r in 1:nperm) {
    for (i in 1:nroi) {
      x.perm[i,r] = x[perm.id[i,r]] #spinning x, spatially permuted bins
      y.perm[i,r] = y[perm.id[i,r]] #spinning y, spatially permuted t-values
    }
  }

## compute mean t-value for each S-A axis bin using spatially permuted y inputs 
### set up spun (y) and empirical (x) df
y.spatialperm <- as.data.frame(y.perm) %>% set_names(sprintf("perm%s", seq(from = 1, to = ncol(y.perm))))
y.spatialperm.empiricalx <- cbind(x, y.spatialperm) %>% as.data.frame()
colnames(y.spatialperm.empiricalx)[1] <- c("empirical.x")

### compute mean t-values for S-A axis bins for all spins
envSES_effects_SAaxisbins_permutedy <- y.spatialperm.empiricalx %>% group_by(empirical.x) %>% #group by empirical S-A axis bins
                                       dplyr::summarize(across(everything(), \(x) mean(x, na.rm = TRUE))) %>% #calculate mean t-value for permuted data
                                       select(-empirical.x) %>% t() %>% as.data.frame() %>%
                                       set_names(c("SA.bin1.tvalue","SA.bin2.tvalue", "SA.bin3.tvalue", "SA.bin4.tvalue"))
enveffects_SAbin_enrichment <- function(bin, enrichment_type){
  
  nperm = 10000
  
  SAbin.tvalue.empirical <- envSES_effects_SAaxisbins %>% filter(SA.axis.quartile == bin) %>% select(mean.tvalue.significant)
  
  #Permutation-based p
    if(enrichment_type == "smaller"){
      perm.p <- sum(envSES_effects_SAaxisbins_permutedy[,bin] < SAbin.tvalue.empirical$mean.tvalue.significant)/nperm
    }
  
  if(enrichment_type == "larger"){
      perm.p <- sum(envSES_effects_SAaxisbins_permutedy[,bin] > SAbin.tvalue.empirical$mean.tvalue.significant)/nperm
  }
  
  return(perm.p)
}

S-A quartile 4

enveffects_SAbin_enrichment(bin = 4, enrichment_type = "larger")
## [1] 0.0323